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ABSTRACT 

We perform two-dimensional hydrodynamic simulations for the thermonuclear explosion of 
Chandrasekhar-mass white dwarfs with dark matter (DM) cores in Newtonian gravity. We include 
a 19-isotope nuclear reaction network and make use of the pure turbulent deflagration model as the 
explosion mechanism in our simulations. Our numerical results show that the general properties of the 
explosion depend quite sensitively on the mass of the DM core Mum- a larger Mum generally leads 
to a weaker explosion and a lower mass of synthesized iron-peaked elements. In particular, the total 
mass of ^®Ni produced can drop from about 0.3 to O.OSMq as Mum increases from 0.01 to O.O3M0. 

We have also constructed the bolometric light curves obtained from our simulations and found that 
our results match well with the observational data of sub-luminous Type-la supernovae. 


1. INTRODUCTION 
1.1. Type-la supernovae 

Type-la supernovae (SNIa) are important astrophysi- 
cal objects b ecause of the similarity in t heir light curves 
and spectra (iBranch &: TammannlflQQ'^ . which leads to 
wide applications of SNIa in cosmological distance mea¬ 
surement such as the determination of the Hubble pa¬ 
rameters (jLeibundgut fc Pintol 119921) and t he discovery 
of the accelerating expansion of the universe (IRiess et al.l 
119981 IPerlmutter et al.lll99'^ . However, despite their im¬ 
portant roles in modern cosmology, both the progeni¬ 
tor system and explosion mechanism of SNIa are not 
yet fully understood. While it is generally believed 
that SNIa are due to the thermonuclear explosion of 
a carbon-oxygen white dwarf (WD) in binary systems, 
it is still unclear whether the companion is a normal 
non-degenerate star or another WD. Traditionally, SNIa 
is attributed to the e xplosion of a WD at the Chan¬ 
drasekhar mass limit (lArnetd I1969D . The WD has a 
mass initially far from the mass limit. Depending on 
the accretion rate, the mass can either gradu ally grow 
until t he baryonic matter becomes degenerate (iNomotol 
I1982all . or a detonation fro nt forms at the envelope and 
sheds away the outer mass (|Nomotolll982bD . Both mech- 
anis ms provide condition s for the formation of a first trig¬ 
ger (|Nomoto et al.l ll 984l) which spread s in the form of a 
deflagration wave (|Nomoto et al.l]T973 l a nd unbinds the 
star. However, neither pure deflagr ation (iNomoto et HI 
I1984D nor pure detonation model (|Arnettl I1969D is ad^ 
equate to explain the observed velocity profile, optical 
light curve, spectra, galactic chemical abundance and ex¬ 
plosion strength. Furthermore, recent studies show that 
SNIa can b e formed without invoking Chandrasekhar 
mass WD (jScalzo et al.l l2014fl . For example, violent 
white d warf mergers can al so explain the SNIa distri¬ 
butions (|Pakmor et al.l[20T^ . 

The difficulties encountered by the pure deflagra¬ 
tion and pure detonation models have led to exten¬ 
sions of models incl uding the pure turbul ent deflagra - 
tion (PTD) model (jReinecke et al.l [T999a]lbl . l2 002albD, 
delay e d-detonatiqn tran sitio n (DDT) model (iKhqholovI 
Il989t iKhokhlo^ ll991allbl H iKhokhlov et al.l I1997D and 


gravitationally confined detonation (CCD) model (previ- 
qusly k nown as the detonation failed deflagration model) 
m ew ^[2?in^ Ik asen fc Plewal 120071: [Jordan et al.l 120081 : 

iMeakin et akll^OOl : l.lordan et al.ll2012fl . Each model has 

its own theoretical difficulties. For example, while the 
PTD mod el can produce exp losion with a variety of 
strengths (iRoeoke et al.l 120061) , there are still unburnt 
low-velocity c arbon and oxygen n ear the core, which are 
not observed (|RoeDke et al'l 20071 1. The DDT model can 
provide sufficient inte rmediate mass elements (IM E) and 
leave very little fuel (|Gamezo et al.l 12004 I2005D . How- 
ever, the possibility of transition is still being debated 
([Lisewski et al.l[2000D . 

In recent years, the PTD, DDT and CCD models 
are studied extens i vely in multi-dimensional simulations 
(iLong et al.l I20l4 iSeitenzahl et al.l 120131 : [Jordan et al.l 
I2012D . The models can well explain the phenomena of 
normal SNIa, i.e., supernovae with a correlated peak lu¬ 
minosity against B-band decline ra te, chemical stratif i- 
cation and a large velocity gradient ([Benetti et al.ll200'^ . 
However, there is a significant number of pecul iar SNIa 
which are sub-luminous and super-luminous (|Li et al.l 
I2001D . In particular, sub-luminous SNIa have a much 
lower absolute magnitude of B- band at maximum. For 
example, the famous SN1991bg (|FiliDDenko et al.lll992D 
recorded a 2.5 mag and 1.6 mag dimmer in B- and 
V-band peak magnitudes. The B-band decline rate is 
faster than the norm, with a lo wer expansion velocit y 
and stronger Si H absorption lines (iDoull fc Baronl2011[l . 
It was initially assumed that such unusual SNIa are ex¬ 
tremely rare. However, there are now sufficient num¬ 
ber of sub-luminous SNIa that they are classified as the 
FAINT group as suggested in ([Benetti et al.ll200^ . De¬ 
tailed study shows that the light curves in this group of 
SNIa are also homogeneous among t hemselves as those of 
normal SNIa (|Doull &: Baronl [20111). The fa intest SNIa 
ever found to date is SN2008ha ( Folev et al.ll2009l) . with 
a low magnitude of My = —14.2 mag and extremely low 
expansion velocity ~ 2000 km s“^. 

In view of the discovery of sub-luminous and super- 
luminous SNIa, the explosion of a Chandrasekhar-mass 
WD can no longer be the sole explanation of SNIa 
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because of the lack of variety in its explosion. The 
sub-Chandrasekhar mass double detonation model is of- 
ten regarded as the expla nation for sub-luminous SNIa 
(jWooslev fc Weav^ I19941 1. The model suggests that 
when the mass accretion of a WD from its companion 
is adequately fast, the matter on its envelope, mostly 
helium, can be ignited and an implosion is triggered 
(|Nomotolll982bf) . The front converges at the WD core 
and a second explosion is created. This model al¬ 
lows a WD to be burnt if the matter is not yet de¬ 
generate. By tuning the host WD mass, less lumi¬ 
nous SNIa can be modeled, which can be fitted to ex- 
plain certain sub-luminous SNIa, for instance SNI99Ibg 
(|Ruiz-LaDuente et al.l 119931) . The helium de tonation is 
found robust in inducing a second explosion (|Fink et al.l 
l2009fl and the pr edicted optical signal is compatibl e 
with observations (jKromer et al.l[21)1(1 : iSim et al.ll2n'T^ . 
However, recent studies of this model with less mas¬ 
sive heliim shell show similar f eature s as normal SNIa 
l|Siml l201d : iRuiter et al.l I2011L 12014) instead of sub- 
luminous ones. Fur thermore, the detonati on might not 
be started robustly (|Livne fc Glasnerlll99r)D . Even when 
a helium detonation is triggered, the detonation wave 
might not penetrate de ep into the carbon/oxygen core 
(|Moll fc WooslevI I2013D , and the distribution of outer 
chemi cal elements can be in conflict with observati on 
data (jEoeflich fc Khokhlovlll996t iHoeflich et al.l[l996D . 

Another popular proposal for explaining sub-luminous 
SNIa is the pure turbulent deflagration model with 
remnant. This model assumes that the deflagration 
only partially burns the WD, and parts of the WD re¬ 
main bounded after the explosion. It is applied to the 
SN2002cx class of the sub-luminous SNIa. The synthetic 
color light curves and the spectra can match well with 
the ob s ervational data (Ijordan et a l.ll20I2tl Kromer et al.l 
l2013al iFink et mi20I4b In (jKromer et ah 120151 ) this 
model is further shown to be in good agreement with the 
faintest SN2008ha. The recent observa tional hints of the 
SN2008ha remnant (|Folev et al.ll20l4) also support this 
model as the origin of this class of sub-luminous SNIa. 

Violent merging of two low-mass WD’s is also a pos¬ 
sible candidate of explaining the sub-luminous SNIa. 
For example, (iPakmor et al.l 120101) showed that this 
model provides a good match to the observe d features 
of SN1991bg, while in (|Kromer et al.l [^013bD the light 
curves and spectra of SN2010Ip are well reproduced. 


1.2. Dark matter astrophysics 

The effects of various dark matter (DM) candidates 
on stellar evolution and structure have been studied 
in details. For example, the effects of DM annihila¬ 
tion as the energy source in ea r ly stars have been con¬ 
sidered in (iSpoiv ar et al. l 120081 Eipamonti et al.l 120101 
Fairbairn et al.l 120081 Fteese et al.l 120091 iSoolvar et al.l 

20091 iHirano et al.l 1201 11) . ~ The DM particl e captu re 
and evaporation rates of th e sun (|Gouldl Il987al lbf) 
and the Earth (|Gou 3 [T^ were studied in early 
1990s. The dense core of compact stars is also a good 
probe of DM (iBer tone fc Eairbairnll2008l : iFan et al.ll2011l : 

Ide Lavallaz & EairbairnI l2010[) . DM particles can an¬ 
nihilate or decay inside a compact star and thus pro¬ 
vide an energy source (iGonzalez fc Reiseneggerl 120101 
iPerez-Garcia fc SilkI l2014b Eor example, it has been 
suggested that the energy is sufficient to maintain 


the surface temperature of WDs (|Moskalenko fc 


2007 : iHooper et al.ll2012D and neu tron stars (|Kouvarii 
2008 ; iKouvaris fc Tinvakovl I2010D . On the other 


hand, non-self-annihilating DM can affect the star 
by its gravity. The self-gravitating DM core inside 
a compact star might coll apse, which forms a black 
hole and engulfs the star (iGoldman fc NussinovI Il989t 
ide Lavallaz fc EairbairnI I2010D . The detection of an¬ 
cient compact stars can thus provide limits on the 
DM scatt ering cross section for differ e nt types of DM 
particles (IKouvaris fc Tinvakovl I2011I iKo uwisI 120121 : 


iMcDermott et al.ll2012l : iBramante et al.l 20131) . 


1.3. Motivation 

Previously, we have studied the equilibrium structure 
and stability of compact stars with cores composed of 
non-self-annihilating DM particles which are modele d by 
an ideal Eermi gas (iLeung et al.l I2011L12012112013D . In 
particular, for DM particle mass of about 1 GeV, we 
found that the DM core can affect the structure of a WD 
significantly. The DM core can be as massive as about 
O.OIMq and the Ghandrasekhar mass limits of these WD 
can be smaller than those without DM by as much as 
40%. An implication of our findings in (jLeung et al.l 
[Ml is that the initial conditions of SNIa might not 
be as universal as generally assumed. In this paper, we 
study how DM affects SNIa explosions by performing 
two-dimensional hydrodynamic simulations of the ther¬ 
monuclear explosions of Ghandrasekhar-mass WDs with 
DM cores. We find that these objects generally have 
weaker explosions and lower masses of synthesized iron- 
peaked elements, and hence they may account for the 
sub-luminous class of SNIa. 

The plan of this paper is as follows: In Section [2] we 
outline the equations and methods that we used in the 
numerical simulations. Section [3] presents the general 
results of our SNIa simulations in terms of the explo¬ 
sion energy, nucleosynthesis, and features of the propa¬ 
gating flame surface. We also compare the bolometric 
light curves constructed from our simulations with the 
observational data of sub-luminous SNIa. Einally, we 
summarize in section |4l 

2. METHODS 

We have developed a two-dimensional hydrodynam- 
ical code with Newtonian gravity to model SNIa. 
The code makes use of the Weighted Essential Non- 
Osdllatory (WENO) sche me for spatial discretization 
([Barth fc Deconinckll 19991) . This is a fifth-order scheme 
which processes piecewise smooth functions with dis¬ 
continuities in order to simulate the flux across grid 
cells with high precision, while avoiding spurious oscil¬ 
lations around the shock. The discretization in time 
is performed by using the five-stage, third-order, non- 
strong stability preserving explicit Runge-Kutta scheme 
([Barth fc Deconin^ I1999D . Various consistency and 
convergence tests have been done to validate our code 
([Leung et al.ll2ni5l) . Here we only outline the essence of 
the code. 


2.1. Initial Model 

In the simulation there are both baryonic normal mat¬ 
ter (NM) and DM. The initial density profiles are ob¬ 
tained by solving the hydrostatic equilibrium equations 
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for both NM and DM: 


dpNM 

dr 

^Pdm 

dr 


G'(Mc(nm) {f ) + -^c(DM) {f )) , . 

--pnm, (1) 

G'(Mc(nm)(’’) + -^c(DM)(^)) , . 

- 2 -PDM, 


where pnm and pdm are the NM and DM density, re¬ 
spectively. The enclosed masses of Mc(nm) and Mc(dm) 
are determined by 


dr 

dMc(BM) 

dr 


= 47rr^pNM, 
= 47rr^pDM- 


(3) 

(4) 


The initial NM is assumed to be isothermal with a tem¬ 
perature of 10® K. The chemical composition is 50% 
and 50% by mass. To construct the initial 
WD models and simulate the NM dynamics, we employ 
the equation of state (EOS) developed and calibra ted in 
(jTimmes fc ArnetUllOOOl : iTimmes fc Swestvl 1 1999( 1. The 
EOS describes the equilibrium thermodynamics proper¬ 
ties of a gas which includes 1 . electrons in the form of an 
ideal gas with arbitrarily degenerate and relativistic lev¬ 
els, 2. ions in the form of a classical ideal gas, 3. photons 
described by the Planck distribution, 4. contributions 
from electron-positron pairs. On the other hand, the DM 
is modeled by an ideal degener ate Fermi gas with a par¬ 
ticle mass moM of 1 - 10 GeV (|Leung et al.ll 201 lL I 2012 L 
1201311 , which is motivated by recent hints on possible de¬ 
tection of GeV scale DM particles in s ome direct DM 
searches, such as DAM A and GoGeNT (|Bernabei et al.l 
120131: lAalseth fc othersi 1201111 . One concern is that the 
admixed DM core may alter the stellar evolution path 
during the main-sequence stage, where deviations from 
standard stellar evolution theory have not been observed. 
We show in the Appendix that the range of Mdm con¬ 
sidered does not bring significant effects on the evolution 
tracks, especially of the carbon-oxygen WD progenitors, 
namely 4-7 Mq main-sequence stars. 

We plot in Fig. [T] the density profiles of two of our 
initial models with wdm = 1 GeV. The upper panel 
of Fig. [T] is a standard model without DM (Model 2D- 
PTD-3-0-c3 in Table [T]). The lower panel shows the NM 
and DM profiles of a stellar model with a DM core of 
mass Mdm = O.O 3 M 0 (Model 2D-PTD-3-3-c3 in Ta¬ 
ble [T|). The central DM density of this model is about 
3 orders of magnitude higher than that of the NM. Due 
to its high compactness, the DM core produces a strong 
gravitational field around the core region and leads to a 
cusp-like structure in NM and steep NM density gradient 
in the core. We remark that t he density profiles shown 
here are different from those in (|Leung et al.ll201^ owing 
to the choice of DM particl e mass and the tot al mass of 
admixed DM. In Fig. 6 of (iLeung et al.ll2013ll the EOS 
of DM particle is chosen to be 10 GeV ideal degener¬ 
ate Fermi gas inst ead of 1 GeV in th is article. Also, in 
the same figure of (|Leung et al.ll201,3h . the admixed DM 
mass is about 1 O“®M 0 , which is one order of magnitude 
lower than those in Fig. [TJ 


2 . 2 . Hydrodynamics 

The thermonuclear explosion of a WD with a DM core 
is inherently a two-fluid system where the NM and DM 
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Figure 1. Upper panel: initial NM density profile of Model 2D- 
PTD-3-0-C3 listed in Table [T] Lower panel: initial NM (dashed 
line) and DM (solid line) density profiles of Model 2D-PTD-3-3-c3. 

couple through gravity. In principle, one has to model 
the dynamics of the two fluids consistently by solving two 
different sets of hydrodynamics equations. However, the 
typical density of DM in our simulations is about two or 
three orders of magnitude higher than that of NM. The 
size of the DM core is also much smaller than the stellar 
radius. As a result, the dynamical time and length scales 
of NM and DM differ by orders of magnitude, and hence 
performing a consistent and accurate two-fluid hydrody¬ 
namics simulation for SNIa would be a computationally 
challenging task. On the other hand, due to its high 
compactness, the dynamics of the DM core is governed 
mainly by its self-gravity. The motion of NM near the 
core is influenced by the DM, but not vice versa. Fur¬ 
thermore, the total energy release and nucleosynthesis 
in the explosion depend mainly on the propagation of 
the flame, which lies well outside the DM core. It may 
thus be reasonable to neglect the motion of DM in the 
explosion. 

As a first step towards understanding the effects of DM 
on SNIa explosions, we only model the dynamics of NM 
in the simulations. The DM core is assumed to be sta¬ 
tionary and affects the NM only through its gravitational 
field. The hydrodynamic code solves the two-dimensional 
Euler equations for NM in cylindrical coordinates (r, z) 
with a detailed nuclear reaction network coupled with 
sub-grid turbulence. The equations are 

|9PNM _ , „ 

^ + v.(H = o, 

(5) 

at 

( 6 ) 

-f p)] — pV ■ V$ -[- (5nuc Qturb 

(7) 

+ V • (pNMqy) = Qturb + V • (pNMl^turbVq), 

( 8 ) 

where pnm, Vr,Vz, Pnm and r are the mass density, veloc¬ 
ities in the r and z directions, pressure and total energy 
density of the baryonic matter. The total energy den- 


— -f V 

dt 


d (pnm?) 
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sity includes both the thermal and kinetic contributions 
T = Pnmc + where e is the specific internal en¬ 

ergy. The specific turbulence energy q is determined by 
Eq. Q . The gravitational potential $ is sourced by both 
fluids and is determined by the Poisson equation 

V^<i> = 47rG'(pNM + Pdm)- (9) 

In Eqs. o and ([5]), Qnuc and Qturb are the heat sources 
from nuclear fusions and sub-grid turbulence, respec¬ 
tively; Qi, is the heat loss due to neutrino emission, and 
is the effective eddy viscosity. We refer the reader 
to (jNiemever fc Hillebrandtl[T995l : lReinecke et iini2002atl 

for a detailed discussion on how these quantities are de¬ 
termined in the simulations. To calculate the heat pro¬ 
duction from nuclear reactions, we incorporate the 19- 
iso tope nuclear r eaction network subroutine developed 
bv lTimmesI (jl999fl into our hydrodynamic code. The iso¬ 
topes include ^H, ^He, ^He, 

28Si, 32s, 36Ar, 40ca, ^SCr, 62Fe, 56^^^ neu¬ 

tron and proton. The fusion network includes reactions 
starting from hydrogen burning up to silicon burning. 
Reactions of (a, 7 ) and (a,p)(p, 7 ) are also included. 

We employ the PTD model as the explosion mechanism 
using the standard configurations that have been consid¬ 
ered in the literature. In particular, an initial flame of 
shape c3 is imposed for all the simulation models and 
the propagation of t he flame is modeled by the standard 
level-set method Isee lReinecke et ahl (jl999a[ l for details). 
Finally, we also construct the theoretical light curves 
from our simulation da ta by using the analytical model 
for SNIa (|Arnettl [198211 . which assumes the photon dif¬ 
fusion limit. This model takes three input parameters: 
the ejecta mass Mej, ejecta velocity Uej and the nickel 
mass Mni, which can be derived from the simulation re¬ 
sults. We also employ the opacity k = O.I g~3c m^ and 
the ga mma-ray deposition function according to lArnettI 

dual. 


3. RESULTS 
3.1. General properties 

In our simulations, the central density of NM is fixed 
to be /Oc(NM) = 3 X 10® g cm“3 because it is expected 
that the minimum density needed for triggering the 
thermonuclear explosion is about 2 — 5 x 10 ^ g cm~3 

(llwamoto et al.lll999l:IWooslevlll997l : lLeSaffre et al.ll^OObt 

iSeitenzahl et al.l 1201 li b We treat the DM core mass 
Mdm as a parameter in the simulations, and we assume 
wdm = 1 GeV unless otherwise noted. 

The properties of four of our typical simulation mod¬ 
els are listed in Table [1] where Pc(NM) (Pc(dm)) and 
Mnm (AIdm) are the central density and total mass of 
NM (DM), respectively. R is the initial stellar radius. 
The first model 2D-PTD-3-0-c3 in the table represents 
a standard model without DM. The other three mod¬ 
els in the table have different DM core masses ranging 
from Mdm = 0.01 to O.O 3 M 0 . All these configurations 
are very close to the corresponding Chandrasekhar-mass 
limits for the given Mdm- The resulting mass of 66 ^ 1 ^ 
energy released through nuclear reactions Enuc, and to¬ 
tal energy Atot are also presented in the table. For the 
same central density of NM, the baryonic mass of the 
star decreases as Mdm increases. However, the (bary¬ 
onic) radius of the star increases with Mdm- 



Figure 2. Upper panel: total energy against time for the models 
listed in Tabled Lower panel: same as above, but for the total 
turbulence energy. 


In the upper panel of Fig.[^we plot the total energy for 
the models listed in Table[I] Note that for a fixed central 
NM density, the total NM mass Mnm decreases as Mdm 
increases. As a result, the initial total energy increases 
with Mdm because most of the binding energy is con¬ 
tributed by NM. At early time, models with less DM have 
faster energy growth than those with a more massive DM 
core. The total energy released, by comparing the initial 
and final energies, decreases when Mdm increases. For 
Model 2D-PTD-3-3-c3, the effects of the admixed DM 
core are so large that the WD remains bound at the end 
of the simulation due to its much lower energy release. 
In the lower panel of Fig. [2l we plot the total turbulence 
kinetic energy against time for the same models. Simi¬ 
lar to the total energy released, the sub-grid turbulence 
energy drops when Mdm increases. The upper panel of 
Fig. [5] also shows that the rate of energy release, reflected 
by the slopes of the curves, decreases as Mdm increases, 
which implies that the time needed for a WD to reach the 
same amount of burnt matter increases. We list the mass 
fractions of major elements at the end of the simulations 
in Table [D In general, the amounts of unburnt fuel and 
IMF increase with Mdmj while those of iron-peaked el¬ 
ements drop. For example, the unburnt are about 
34% and 46% of the total mass for models 2D-PTD-3-0- 
c3 (Mdm = 0) and 2D-PTD-3-3-c3 (Mdm = O.O 3 M 0 ), 
respectively. On the other hand, the mass fraction of 
66Ni decreases significantly from about 24% to 2.8% as 
Mdm increases from 0 to O.O3M0. This is related to the 
different initial density distributions in the models. For 
a larger Mdm, the amount of matter that can reach suf¬ 
ficiently high density for complete combustion decreases. 

Next we consider the effects of DM on the flame sur¬ 
face. We plot in Figs. [5115] the flame surfaces (represented 
by the temperature) at f = 1 s for the models listed in 
TablelD In Model 2D-PTD-3-0-c3, which is the standard 
PTD model without DM considered in the literature (see 
for exa mple INiemever &: Wooslevl (1199711 , IReinecke et al.1 
(|1999all and IReinecke et al.l ( 2002alP . the flame shows 
a convoluted structure, with clear instabilities of flame- 
fluid interaction including the Rayleigh-Taylor instabili¬ 
ties and Kelvin-Helmholtz instabilities. The injection of 
fuel into the flame can be seen as well. When Mdm in¬ 
creases, the injection of fuel can still be found. But the 
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Table 1 

Simulation setup for four PTD models: central densities of NM Pc(NM) DM Pc(DM) s-''® in units of 10® g cm~®. Masses of the baryonic 
matter MnM) dark matter Mdmi nnd the final nickel-56 mass Mn; are in units of solar mass. R is the initial stellar radius. E^uc and Etot 
are the energy released by nuclear reactions and final total energy, respectively, both in units of 10®® erg. We assume rriDM = 1 GeV. 


Model 

Pc(NM) 

Pc(DM) 

Mnm 

Mum 

R (km) 

Mni 

-F'nuc 

Etot 

2D-PTD-3-0-C3 

3.0 

0.0 

1.377 

0.00 

1.92 X 10® 

0.33 

7.2 

2.1 

2D-PTD-3-1-C3 

3.0 

150 

1.313 

0.01 

2.22 X 10® 

0.25 

5.8 

1.5 

2D-PTD-3-2-C3 

3.0 

530 

1.223 

0.02 

2.79 X 10® 

0.14 

3.9 

0.46 

2D-PTD-3-3-C3 

3.0 

1050 

1.015 

0.03 

4.25 X 10® 

2.9 X 10-® 

1.1 

-0.85 


Table 2 

Mass fractions (normalized by the stellar mass) of all isotopes at the end of the simulations for the models listed in Table [T] 


Isotope 

2D-PTD-3-0-C3 

2D-PTD-3-1-C3 

2D-PTD-3-2-C3 

2D-PTD-3-3-C3 

I^C 

0.30 

0.32 

0.36 

0.43 

16 o 

0.33 

0.35 

0.39 

0.46 

24Mg 

1.4 X 10-2 

1.4 X 10-2 

1.4 X 10-2 

1.0 X 10-® 

28Si 

5.2 X 10-2 

5.9 X 10-2 

4.9 X 10-® 

3.9 X 10-2 

®2s 

1.8 X 10-2 

1.9 X 10-2 

1.6 X 10-2 

1.1 X 10-2 

®®Ar 

3.2 X 10-® 

3.1 X 10-® 

2.5 X 10-® 

1.7 X 10-® 

«Ca 

3.1 X 10-® 

2.9 X 10-® 

2.3 X 10-® 

1.3 X 10-® 

44 Ti 

6.5 X 10-5 

9.6 X 10-® 

7.3 X 10-® 

2.4 X 10-® 

48 Cr 

1.5 X 10-4 

1.4 X 10-4 

8.8 X 10-4 

2.3 X 10-5 

52 Fe 

4.8 X 10-® 

3.8 X 10-® 

2.3 X 10-® 

5.1 X 10-4 

54 Fe 

5.1 X 10-2 

4.5 X 10-2 

3.3 X 10-2 

1.7 X 10-2 

56Ni 

0.22 

0.18 

0.11 

2.8 X 10-2 


Kelvin-Helmholtz instabilities are suppressed. The flame 
surface also becomes smoother and less turbulent. 

3.2. Connection with Sub-luminous Supernovae 

The light curve of a sub-luminous SNIa has a low 
peak luminosity, suggesting that the ®®Ni content is 
lower than ordinary SNIa. The explosion is very weak 
and only partial ejecta are dispelled instead of a dis¬ 
ruption of the whole star. Notice that, how much and 
how the mass is ejected in the explosion are not clear 
unless the simulation continues until the homologous 
expansion phase is reached, which takes place about 
ten sec onds after the deflagration/detonation st age has 
ended (iRoeoke fc HillebrandtJl20n^ lRoeDkell2005H . How¬ 
ever, this involves using either a sufficiently large simula¬ 
tion box which can accommodate the rapidly expanding 
ejecta or expan ding meshes t o prevent the ejecta from 
leaving the box (jRoeDkell2005[) . 

Since the amount of ejecta mass is an important pa¬ 
rameter in constructing the resultant light curves from 
the simulations, we constrain it in the following ways. 
First the minimum ejecta mass is estimated by counting 
all the fluid elements with positive energy at the end of 
the simulation. Second we assume that the maximum 
ejecta mass is equal to the total mass of the star. 

In Fig. Owe show the bolometric light curves for the 
four simulation models listed in Table [T] (from the top 
solid line to the next-to-bottom solid line) and also one 
additional model (the bottom line) not listed in the table. 
The extra model has a DM core mass Mdm = O.O32M0. 
For the first three models (2D-PTD-3-0-c3, 2D-PTD-3-1- 
c3, 2D-PTD-3-2-c3), the total final energies are positive 
and hence we use the maximum ejecta mass for each 
model to construct the light curves. On the other hand, 
we use the minimum ejecta mass to construct the light 
curves for the remaining two models because their total 
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Figure 3. The flame surface of Model 2D-PTD-3-0-c3 at £ = 1 s. 
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Figure 4. Same as Fig. [3l but for Model 2D-PTD-3-l-c3. 
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Figure 5. Same as Fig. [S] but for Model 2D-PTD-3-2-c3. 
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Figure 6 . Same as Fig. O but for Model 2D-PTD-3-3-c3. 

final energies are negative. Fig. 0 shows that the peak 
luminosity depends sensitively on Mdm- In particular, it 
can decrease by almost two orders of magnitude as Mdm 
increases from 0 to O.O32M0. In the figure, we also plot 
the data from the constructed bolometric light curves 
for some examples of sub-luminous SNIa for comparison. 
It can be seen that our SNIa simulations with admixed 
DM give a range of peak luminosities that covers the 
observed sub-luminous SNIa including the exceptionally 
dim SN2008ha. 

Our results suggest that the variations of the observed 
light curves of different sub-luminous SNIa may be due 
to the fact that the underlying WDs of the systems con¬ 
tain different amounts of DM. For a given observed SNIa 
light curve, we can use Mdm as a parameter for perform¬ 
ing hydrodynamical simulations to fit the observed data. 
In principle, for a given Mdm, a unique light curve can 
be determined by the resulting velocity profile and total 
mass of ejecta obtained from the simulation. However, 
as we discussed above, the total mass of ejecta cannot be 
determined accurately from the simulations due to com¬ 
putational limitations. We thus calculate two different 
light curves corresponding to the minimum and maxi¬ 
mum ejecta masses for a chosen Mdm- 

In Fig. [5] we plot the bolometric light curves by us¬ 
ing the results from a simulation model with Mdm = 
O.O2 M0. The solid and dashed lines in the figure are 



Figure 7. From the top solid line to the next-to-bottom solid 
line: bolometric light curves of Models 2D-PTD-3-0-c3, 2D-PTD- 
3-l-c3, 2D-PTD-3-2-C3 and 2D-PTD-3-3-c3. The lowest light curve 
corresponds to an extra model, with similar configurations as the 
above four models but with Mqm = 0.032Mq. Observational 
data of SN2011ay, SN2005hk, SN1999by, SN2003gq, SN2005cc and 
SN2008ha are included. 



Figure 8. Bolometric light curves of a simulation model with 
Mdm = 0.02Mq. The solid (dashed) line is constructed by assum¬ 
ing that the ejecta mass takes the minimum (maximum) value esti¬ 
mated from the simulation as discussed in the text. Observational 
data of SN1999by is also included. The error bars correspond to 
the uncertainties in the distance modulus and the measurements. 

obtained by using the minimum and maximum ejecta 
masses, respectively. This simulation model has a pos¬ 
itive total final energy. As a result, the minimum and 
maximum ejecta masses estimated from the simulation 
are comparable, and hence the two constructed light 
curves are also quite close to each other. We expect that 
the light curve corresponding to the actual ejecta mass 
should lie between the two limits. In the figure, the ob¬ 
servational data of a sub-luminous supernova SN1999by 
(with error bars) are also presented for comparison. The 
error bars correspond to the uncertainties in the distance 
modulus and measurements. It is seen that the data 
around the peak luminosity lie very close to the region be¬ 
tween the two constructed theoretical light curves, which 
represents effectively our uncertainty in the calculation. 
At later time, the observational data decays faster than 
the theoretical light curves. One possible reason may be 
that our assumption of the opacity law and gamma-ray 
deposition function are no longer valid at later time. 
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Figure 9. Same as Fig. [8] but for a simulation model with 
-^DM = 0.026AfQ. Observational data of SN1991bg is included 
for comparison. The error bars correspond to the uncertainties in 
the distance modulus and the measurements. 

As a different example, we plot in Fig. the bolo- 
metric light curves of a simulation model with Mdm = 
O.O26M0. Contrary to the simulation model with 
ATdm = O.O 2 M 0 discussed above, this model has a neg¬ 
ative total final energy and hence the estimated mini¬ 
mum ejecta mass differs from the maximum ejecta mass 
quite significantly. As a result, the two corresponding 
light curves (solid and dashed lines) are not close to each 
other. In the figure, the observational data of another 
sub-luminous supernova SN1991bg is also plotted, which 
agrees quite well with the theoretical light curve con¬ 
structed with the minimum ejecta mass. 

The above examples show that the admixture of dark 
matter can produce SNIa light curves with large varia¬ 
tions in peak luminosities comparable with those of sub- 
luminous SNIa. Mej also affects the peak luminosity, but 
its influence is much less pronounced than Mdm- On the 
other hand, Mej dominates the width of a light curve. 
This suggests that the observational data of an SNIa can 
provide hint on the ejecta mass and its admixed DM 
mass, with Mdm determining the peak luminosity while 
Mej the light curve width. Given the light curve data 
of an SNIa, we search for the best-fitted pair of Mdm 
and Mej, where the values of Mni and z;ej are derived 
from simulations as functions of Mdm- Models with the 
minimum chi-squared values are chosen to be the rep¬ 
resenting models of that SNIa. In Table [3] we list the 
models with minimum chi-squared values of several well 
observed sub-luminous SNIa. The table lists the agree¬ 
ing Mej, Mdm with their implied Mni, ejecta velocity Uej 
and the two ejecta mass limits Mgj (max) and Mgj (min); 
which are derived from simulations. We regard that the 
DM admixture of this model can be a possible explana¬ 
tion of an observed sub-luminous SNIa if its fitted Mej 
lies within the two limits of Mej. All sub-luminous SNIa 
in the list except SNI99Ibg and SN1993H give an Mej 
inside the two limits, showing that these SNIa could pos¬ 
sibly have admixed DM cores in the WD progenitors. 
The SN1993H (SN1991bg) has an ejecta mass just above 
(below) the upper (lower) Mej limit, showing that the 
observational data is declining slower (faster) than what 
the current configuration can provide. 

In Fig. [To] we plot the peak-luminosity against the fit¬ 
ted Mdm of the mentioned SNIa. The peak luminosities 


Figure 10. The SNIa peak-luminosity against Mdm- Peak lu¬ 
minosities of some observed SNIa are also plotted. The shaded 
regions are excluded by the mass bounds of the progenitor mass 
and minimum ejecta mass derived from simulations. 



time (s) 


Figure 11. Total energy against time for Models 2D-PTD-3-0-c3, 
2D-PTD-3-1-C3 and an extra model similar to Model 2D-PTD-3- 
l-c3 but with Mdm = W~^Mq and uidm = 10 GeV. 

of some observed SNIa are included as data points. The 
shaded regions of the plot are excluded by this model be¬ 
cause these regions correspond to models with an ejecta 
mass out of the bounds. Given an Mdm the range of 
peak luminosities is very small compared to the observed 
range of SNIa peak luminosities. This implies that the 
admixed DM mass can be well constrained by the peak 
luminosity. 

In summary, our work shows that the admixture of 
DM can explain the observation data of sub-luminous 
SNIa. However, as discussed in Sec. I, it should be 
noted that matching of the bolometric light curves of 
sub- luminous SNIa can also be achieved in ot her models 
(iPakmor et al.l 120101 iKromer et al.l l2013allbl: iFink et al.l 
l2014HKromerLtalMn5lb 

4. CONCLUSION 

In this paper, we have performed two-dimensional 
Newtonian hydrodynamic simulations to study the ef¬ 
fects of DM on the thermonuclear explosion of WDs near 
the Chandrasekhar mass limit. Our initial models are 
constructed by solving the two-fluid hydrostatic equilib¬ 
rium equations for NM and DM with a fixed central NM 
density of 3 x 10® g cm“^, which is expected to be near 

































Table 3 

Fitting results of observed sub-luminous SNIa. Masses are in units of solar mass and the ejecta velocity v^j are in units of 10® cm s“^. 
Mgj (max) S'lid Mgj (min) ^^e the maximum and minimum ejecta masses derived from simulations. The last column marks the possibility 
of using admixed DM to explain the observed SNIa. moM = 1 GeV is assumed. 


Supernova 

Mdm 

Mej 

Mni 

fej 

-^ej (max) 

-^ej (min) 

DM origin 

SN1991bg 

0.025 

0.20 

0.084 

0.44 

1.14 

0.21 

No 

SN1993H 

0.008 

1.37 

0.216 

0.74 

1.35 

0.69 

No 

SN1999by 

0.019 

0.61 

0.148 

0.58 

1.24 

0.54 

Yes 

SN2003gq 

0.024 

0.82 

0.094 

4.61 

1.16 

0.25 

Yes 

SN2005CC 

0.027 

0.33 

0.058 

3.70 

1.01 

0.14 

Yes 

SN2008ae 

0.022 

0.67 

0.120 

5.22 

1.20 

0.38 

Yes 

SN2008ha 

0.032 

0.12 

0.004 

2.19 

1.00 

0.09 

Yes 

SN2011ay 

0.000 

1.35 

0.320 

7.63 

1.38 

0.73 

Yes 


the minimum density for triggering the explosion. The 
typical models studied by us are solar-mass WDs with 
small DM cores (~ O.OIM 0 ) formed by DM with a parti¬ 
cle mass of 1 GeV. As a first step towards understanding 
the effects of DM on SNIa, we assume that the DM core 
is stationary during the evolution, and we only model 
the dynamics of the NM fluid. This should be a good 
approximation as the DM core should be affected mainly 
by its self-gravity due to its high compactness. We em¬ 
ploy the PTD model as the explosion mechanism and use 
the standard level-set method to model the flame surface 
during the dynamical evolution. 

We have only considered the PTD model as the ex¬ 
plosion mechanism in this work. It would be interest¬ 
ing to extend our work by using other possible explosion 
mechanisms such as the DDT and GCD models as dis¬ 
cussed in Sec. [T] Finally, we have also assumed that 
the DM core is formed by non-self-annihilating fermionic 
DM with particle mass 1 GeV. How would our results be 
changed if one considers much more massive DM particle 
candidates? We plot in Fig. [TT]the total energy against 
time for Models 2D-PTD-3-0-c3, 2D-PTD-3-I-c3 and an 
extra model similar to Model 2D-PTD-3-l-c3 but with 
AIdm = 1O“^M0 and ttidm = 10 GeV. The initial mass 
of the extra model is 1.35 Mq. Despite the small Mdm 
in the new case, the effects of DM on the energy release 
during explosion is comparable with Model 2D-PTD-3-1- 
c3, w hich has Mdm = O.O lMf:^. As found in our previous 
work iLeung et al.l (|2013[ 1. for the same Mdm or central 
density, a higher todm has a stronger effect on the den¬ 
sity profile comparing to the case of wdm = 1 GeV. It 
is because the DM core is more compact and creates a 
stronger gravitational attraction field, which changes the 
initial density profile more significantly. These changes 
in the density profiles are also reflected by the drop of 
energy release in SNIa explosions. Our results show that 
SNIa explosions are sensitive to both the drop of Ghan- 
drasekhar mass, as shown in the 1 GeV case, and to 
the particle mass. However, modeling more massive DM 
core with moM = 10 GeV is difficult because the region 
expected to be admixed with DM will be even smaller, 
implying that a much higher resolution is needed in order 
to model both DM and NM consistently. 

Our numerical results show that an increase in Mdm 
or Todm leads to a change in the SNIa explosion by either 
decreasing the Ghandrasekhar limit for low todm or alter¬ 
ing the density profile for high ttidm- First, the explosion 
becomes weaker and the total energy release is reduced. 
The total turbulence energy, which is an important in¬ 


dicator for the PTD model, also decreases. Second, the 
amounts of unburnt fuel and IMF increase, while those 
of iron-peaked elements decrease. In particular, the to¬ 
tal mass of ®®Ni depends quite sensitively on Mdm and 
decreases from about 0.3 to O.O3M0 as Mdm increases 
from O.OI to O.O3M0 for the ttidm = 1 GeV case. Finally, 
the Kelvin-Helmholtz instabilities are suppressed and the 
flame surface also becomes smoother and less turbulent 
as Mdm increases. We have also constructed the bolo- 
metric light curves from our simulations and compared 
them with the observational data of sub-luminous SNIa. 
Our results shows that varying the DM core mass from 
about 0.01 to O.O 3 M 0 yields a range of peak luminosi¬ 
ties that covers the observational data very well. The 
variations of the observed light curves of different sub- 
luminous SNIa may be due to the fact that the precursor 
WDs contain different amounts of DM. 
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6 . APPENDIX: EFFECTS OF DM ADMIXTURE ON SNIA 
PROGENITORS 

In this article we have studied how the gravity of DM 
affects the explosion energetics of SNIa. We have shown 
that with an admixed DM core with a total mass in the 
order of the ^®Ni production can be signifi¬ 

cantly suppressed and the corresponding light curves are 
comparable with those of sub-luminous SNIa. However, 
it remains unclear whether such a DM admixture can 
leave observable consequences already during the main- 
sequence phase, which is well constrained by observa¬ 
tional data. Therefore, it is necessary to check if stars 
with DM admixture have unusual evolution paths that 
are inconsistent with observational data, and if the chem¬ 
ical compositions of the resultant white dwarfs are dif¬ 
ferent from those of conventional cases. 

A star acquires DM mainly by accretion through DM- 
NM scattering or by its inherent admixture that exists 
already during its formation stage, where DM acts as a 
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stellar seed. However, following (jKouvarisI 1200811 to es¬ 
timate the DM accretion rate, using conventional DM 
parameters, the typical DM accretion rate is insignifi¬ 
cant compared with the original mass of the star, even 
when we consider a duration of cosmological timescale. 
Therefore, it is unlikely that a star can acquire DM with 
a mass comparable with the host simply by accretion. 
We thus focus on DM which acts as a stellar seed. In 
that case, the gravity of DM is important even in the 
protostellar phase. 

We performed main-sequence star simulations by using 
an open-source stellar evolution code ME SA (Modules 
for Experiments in Stellar Astrophysics) (jPaxton et al.l 
1201 IL [ 201 ^ 120151 1. which can follow the evolution of a star 
from the protostellar phase up to the white dwarf stage. 
We used the MESA code version 3372, which solves the 
fully coupled one-dimensional structure and composition 
equations simultaneously, using the Helmholtz EOS to 
describe the thermodynamics properties of NM. The DM 


9 








Table 4 

Stellar properties at the end of simulations and main-sequence lifetime for models with Mnm = 4 or 7Mq. Mjje and Mqq are the masses 
of ^He, and in units of solar mass. Xc and Xq are the mass fractions of and of the innermost mass shell. Tjj start 
(7 h end) and Tjje start (THe—end) are the beginning (ending) time for hydrogen amd helium burning in units of years. 


Mnm 

Mdm 

Mhb 

Mco 



Th start 

Th end 

TjJe start 

^He end 

4 

0 

0.24 

0.18 

0.40 

0.58 

1.12 X 10*^ 

1.49 X 10® 

1.55 X 10® 

1.89 X 10® 

4 

0.01 

0.24 

0.19 

0.47 

0.50 

1.08 X 10® 

1.46 X 10® 

1.50 X 10® 

1.82 X 10® 

4 

0.02 

0.18 

0.37 

0.56 

0.42 

1.05 X 10® 

1.44 X 10® 

1.47 X 10® 

1.82 X 10® 

4 

0.03 

N/A 

N/A 

N/A 

N/A 

1.01 X 10® 

1.44 X 10® 

1.47 X 10® 

N/A 

7 

0 

0.32 

0.13 

0.44 

0.54 

2.30 X 10® 

3.98 X 10'^ 
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Figure 12. The H-R diagrams of stars with 4 Mq NM and 0.03 
Mq dm but with different innermost mass shell. The simulation 
is done until the time-step becomes smaller than 10^ years. 

is assumed to be in hydrostatic equilibrium, and to a 
good approximation, the DM profile remains static dur¬ 
ing the simulation. We observed that due to the com¬ 
pactness of the DM core, in most of the stellar lifetime, 
the DM core has a size smaller than the outer radius of 
the innermost fluid elements. Effectively, we modified 
the hydrostatic equation in the MESA code by including 
the DM core which behaves like a point-mass as 

dP^M G(wnm(?') + AfpM) , 

dm 

All notations have the same meaning as those in the main 
text. Due to the singular behavior of the DM poten¬ 
tial near the core, there are numerical difficulties that 
the results are resolution dependent. Also, the typical 
time-step becomes prohibitively small, due to the large 
potential gradient, which leads to a large density gra¬ 
dient and hence a large chemical composition gradient 
near the core. Also, the 1/r potential from the DM leads 
to divergence in constructing the initial model. To ame¬ 
liorate these problems, we smear out the effect of DM 
by increasing the innermost fluid elements from 10“® to 
1 O“^M0. This allows us to capture the effects of DM’s 
gravity within a reasonable simulation time. 

We use the star evolution model lM_pre_ms_wd in the 
test suite package to follow the stellar evolution from the 
protostellar phase. We considered star models with a 
mass from 4 to 7 Mq, which are believed to be the pro¬ 
genitors of carbon-oxygen white dwarfs. 

To show that the rise of innermost fluid element mass 



Figure 13. The H-R diagrams of stars with 4 Mq NM. Mdm 
ranges from 0 to 0.03 Mq. All simulations have the innermost 
mass shell fixed at 10 “^Mq, and they are stopped when the time- 
step becomes smaller than 10^ years. 



Figure 14. Same as Fig. 1131 but for Mnm = 7Mq. 


does not introduce spurious results, we plot in Fig. [T^] 
the HR diagram of a star with 4 Mq NM and 0.03 Mq 
DM, but with different innermost fluid element masses. 
We do not follow the whole evolution till the formation 
of the white dwarf because the timestep has already be¬ 
come prohibitively small when it enters the helium burn¬ 
ing phase. We stopped the simulation when the average 
time step drops below 10^ years. In the mass range con¬ 
sidered, the qualitative behavior of the stellar evolution 
remains unchanged. This shows that in this resolution 
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the basic properties of the niain-sequence phase are cap¬ 
tured. The model with a higher innermost mass shell 
can run longer due to the stronger smearing of the DM 
point-mass gravity. 

We plot in Figs. [13] and |T4| the HR diagrams for 
star models with a normal matter mass of 4 and 7 solar 
masses, but for different Mdm- In both figures, the evolu¬ 
tion paths of the hydrogen burning phase and the helium 
burning phase are insensitive to Mdm- We terminated 
the simulations for Mdm = 0.03M before the exhaus¬ 
tion of core helium because of the small time-steps. One 
qualitative difference that can be observed is the disap¬ 
pearance of the horizontal branch during helium burning 
for the model with Mnm = 7Mq and Mdm = O.O 2 M 0 . 

In Table m we tabulate the stellar properties extracted 
from profiles at the end of simulations, and also the age of 
the star when hydrogen or helium burning starts or ends. 
No results are listed for models with Mdm = O.OSMq 
because the simulations are terminated before the helium 
burning phase commences. For models with Mnm = 
4 M 0 , when Mdm increases, the helium mass decreases 


while the carbon-oxygen mass increases. Also, hydrogen 
burning begins and ends sooner, with the whole hydrogen 
burning lifetime shortened. Similar features are observed 
for the helium burning. The mass fraction increases 
while that of decreases. The effects of DM in models 
with Mnm = 7Mq become smaller so that there is almost 
no change when Mdm = O.OIM 0 . But as Mdm further 
increases, similar effects can still be observed, including a 
lower helium mass, earlier hydrogen and helium burning 
and a shorter main-sequence lifetime. Also, an increase 
(a decrease) in (^® 0 ) mass fraction is observed. 

From the above comparison, we have shown that in 
the mass range of Mdm considered in the main text, the 
DM core which is assumed to exist as early as the star 
forms, does not alter the stellar evolution significantly. 
Specifically, all models predict a path during the main- 
sequence phase in the HR diagram comparable with the 
cases without DM. Moreover, the final chemical compo¬ 
sition of the carbon-oxygen white dwarf does not deviate 
significantly from what we have assumed, a 50 % carbon 
and 50 % oxygen by mass. 
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